{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Predictions with Pyro + GPyTorch (High-Level Interface)\n",
    "\n",
    "## Overview\n",
    "\n",
    "In this example, we will give an overview of the high-level Pyro-GPyTorch integration - designed for predictive models.\n",
    "This will introduce you to the key GPyTorch objects that play with Pyro. Here are the key benefits of the integration:\n",
    "\n",
    "**Pyro provides:**\n",
    "\n",
    "- The engines for performing approximate inference or sampling\n",
    "- The ability to define additional latent variables\n",
    "\n",
    "**GPyTorch provides:**\n",
    "\n",
    "- A library of kernels/means/likelihoods\n",
    "- Mechanisms for efficient GP computations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import pyro\n",
    "import tqdm\n",
    "import gpytorch\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we will be doing simple variational regression to learn a monotonic function. This example is doing the exact same thing as [GPyTorch's native approximate inference](../04_Variational_and_Approximate_GPs/SVGP_Regression_CUDA.ipynb), except we're now using Pyro's variational inference engine.\n",
    "\n",
    "In general - if this was your dataset, you'd be better off using GPyTorch's native exact or approximate GPs.\n",
    "(We're just using a simple example to introduce you to the GPyTorch/Pyro integration)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x11ddf7320>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANwAAACbCAYAAAAJB/VVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAARHElEQVR4nO3df4yU9Z3A8fcDsyysIAsoriwLrphN0bZeyKeh0VIFm56/Em2ldk3OH40x2k3Os5xZ4ErhIBo5juB5NdpEiT32omuNYZVeo7XS1mqM9oP2qMVzdwVZgt2AXRZ/gLvO7nN/zDMwuzvzzM7OPM88M/N5JQTm+8w889mH/czzfH88n3Fc18UYE45JxQ7AmEpiCWdMiCzhjAmRJZwxIbKEMyZElnDGhMgSzpgQxYLcuYh8E9igqleMar8EuJREwj+hqkfStQUZmzHFEGjCqeorIjItzaYHgMuBBmAj8MMMbWmtWbPGZutNZG3evNnJtC3QhPMMpj4QkSYgrqou0CMiy9K1Zdvpxo0b07YfOXKEuXPnFiDswopqXBDd2Eoxrg0bNvi+NoyEG+0s4FjK41kZ2nwdOZL+irO/vz+f2AIT1bggurFFLa6dO2vYsmUmH344n3nzhmhtPc53vnMip30UI+GOAjUpjwcytPny++SL4qciRDcuiG5sUYmrvX0Sa9fGOHEicbV4+HCMtWtnM3PmmTQ3D497P6ElnIhMBmpUtUtEpnpt5wO/S9eWy77j8Ti9vb0MDAwwNDTExx9/XOjw8zY8PFy0uBzHobq6mrq6OmKxYnzGlr71608nW9KJEw7r18dobh7M8Kqxgh6l/AqwSES+DDQC3wBWA/eJSCswFfix9/R0bePS29vLGWecQX19PfF4nKqqqsL9EAXyxRdfFC0u13Xp7++nt7eX+fPnFyWGUnfoUG7tmQQ9SvlnEqOOAO8Au7z23cDuUc8d0zZeAwMD1NfX4zgZB4cqmuM41NbW0tfXV+xQSlZDA/T0pG/PRVlMfLuua8mWheM42L2PE7dpU5yampHHr6bGZdOmeE77KYuEG6/29kk0NU1h2rQpNDVNob099x9/3bp1tLW1cfHFF9PW1kZra6vv8zMNE2/fvp2DBw/m/P6HDx/mlltuyfl1Jj/NzcM88kicBQtcHMdlwQKXRx6J5zRgAsUZpSyK9vZJtLSc7vj29EBLSwzI7aCtWrWK2bNns2PHDm6++WaOHTvm+/xM84W33377uN8zVX19PVOmTBnTrqoMDg5yySWXTGi/Jrvm5mGamwfzmh+smIQr1CjT7NmzRzyeNWsWGzduxHVdampq+OSTT7jwwgs5fvw411xzDevWrWPz5s3cddddXHvttbz00ks8+eSTtLS0sHr1arZt28aSJUvYuXMn27ZtY2hoiBdeeIGuri6qqqp46KGHgEQ/9dFHH6WhoYGPPvqI/v5+tm7dSnV1NVdddRUdHR2ceeaZ1NbW0tHRQV9fH3fffTcLFizI/+CZgqmYS8pCjTKl09jYyLJly2htbeX6669n7ty57Nmzh4aGBqqqqjj33HOZPXs2d9xxB/39/cRiMRoaGpg8eTJ1dXVcccUVrFy5ElVl9+7dLFq0iKuvvprFixefeo9nnnmGpqYmbrjhBubMmcPUqVO58sormT59Ou+++y5NTU0sXbqUuro6LrvsMqqqqvjggw/y/+EqSCG6HNlUTMJlGk3KdZQpk1gshuu67Nixg8svv3zMAEVyUCfTwMWkSZOIx+OsXLmSvXv3MjAwwK233npq+7Fjx06NMjqOw+HDh3n99ddZsmTJiEGj9957j+7ubhYvXmyDJDlIdjl6ehxc16Gnx6GlJVbwpKuYhCvUKBNAV1cXBw4c4P333wegu7ubvXv3ArBv3z7a2to4evQoR48epbu7m4MHD3LgwAF6enro6+s79fquri7279/P/v37T+3v7bff5o033qCjo4OOjo5T73njjTfy9NNPs337dnp7e+nr6+Pll1/mrbfeYs+ePZx33nk899xzxONxdu3aRXd3N2+++WYeR6yy+HU5CskpxU/BNWvWuKmDEd3d3VxwwQWA/wRze/sk1q+PcehQ4sy2aVPuo0wTNd6J78cee4zbbruNWCzG888/z3XXXVewGFKPU6pSXCRcaNOmTcF1x04tOY7LyZMj+/jZFi8X+26ByEiOMkXZZ599xr333stFF13E0qVLix1OxSjUxHY2FZVwpeCee+4pdggVadOm+IhpI5h4l8NPWfThbBVFdrYax1+hJrazKYszXHV1Nf39/dTW1hY7lEhKLl6urq4udiiRFkaXoywSrq6u7tTI3dDQEJMmRe/EPTw8XLS4Um/PqWTFHDRLKouEi8Vip247sRE3k06hlvblK3qnAmMCENY8WzaWcKYiBLm0LxeWcKYiBL20b7ws4UxFKOTSvnwEXdPkn4EjwExVfTil/UXgS4ALuKra6LV3AF8HdqnqHUHGZipLYmAkXvRRysDOcCLyDWCOqrYBs0Rkqdc+A1ilqgtJJN1TXvvXgEdVtc6SzQShuXmYzs5BTp4cpLNzMPRkg2DPcFcD73r/3uc9fkNVPwH+4rV/G/i19+/lwD+KyG7gh6rqW2HTCsEWTlRjK8e4gky41GrKnwPpZl2XAWsAVHWLiGwD/s1rW++3cysEW1hRjS3XuMKa3J7o8Qpy0CS1mvIM4G+pG0UkBgyp6lCyTVXjJOpWNgYYlylTYd1Emo8gI/kV8FXv3xcCL4rIzJTty4HfJh+ISHJWcgbwaoBxmTIVlcltP4ElnKq+BnwuIj8A+r0/P0t5ynJGFn59VUQeBr4LPB5UXKZ8RWVy20/QlZfvG9V0U8q2fxn13EuDjMWUv7BuIs1HdC5ujclTVCa3/VjCmbIR1k2k+bCEMyUlW+3IKExu+4nO8I0xWUTlnrZ82BnOlIxSGPbPxhLOlIxSGPbPxhLOlIyo3NOWD0s4UzJKYdg/G0s4UzJKYdg/m9LpbRpDaZSr92NnOGNCZAlnTIgs4UzkJFeTLFw4P7BvIi0W68OZSCmH1SR+yuejw5SFclhN4scSzkRKOawm8WMJZ0Lnt+K/HFaT+ClKIVhv24iiryLSBHwfOOG1dQYZmymObH20sL6JtFhCLwTrbUtX9PUh4EHgYWBzUHGZ4srWRyuH1SR+Qi8E6z0eUfSVRMnzRar6KYCINIpIzCubZ8rIePpoydUk5fideuNKOMdx/h14wHXdvhz2nbEQbJqirz8DPk55bRw4G/hrpp1b5eXCCTO2efPO5fDhsb928+YNjfk/jeoxC6Py8lPA3zuOczbwf8BLbvZvsfctBKuqcRFZDTzhbZuasrmGRFm9jKzycmGFFdv997u0tLhj+mj33++mjSGqxyzoyst/Bp4FBoD/BO5zHOcux3Fm+LwmYyHY0UVfVXUAOCgiNSIyFTikqidz/FlMCSj3Plo24z3D/QaYRuJMJ67rfuo4zjTgf4AV6V6gqq+JyPI0hWBvIlH09W3gbU4XfV0NtJJI6lUT/HlMCSj1Ff/5GG/C/QH4V9d1Uwcx4sD/+r0oUyHYdEVfVfUd4J1xxmNMSRrXJaXruutGJRuu637huu6PggnLlLJspewqWXksUDORUe6Lj/NlHz2moMp98XG+LOFMQZX74uN8WcKZCcnUTyv3xcf5svO8yZlfP63cFx/nyxLO5Myvn9bZOQjEQ/me7VJkCWdylq2fVskT29lYH87kzPppE2cJZ3JWDiXHi8USzqTlt1qk0hcg58P6cGaM8awWsX7axNgZzoxhq0WCYwlnxrDVIsGxhDNj2ChkcCzhzBg2ChkcSzgzho1CBsd6wSYtG4UMRlEqL4vITcA9wJnAzaqqXvtPge8Bf1LVK4OMzSSG/xNrHufbmseQhF552avYdUJVlwJbgY1eez3wlleN2ZItYMm5tp4eB9d16OlxaGmJWTmEgAV5dNNVXkZVXVV9zmv/I6eLva4AfiIivxSRswKMy2BzbcUS5NHNWHk5xbeAbQCq2iYi/03iUnMbcIvfzq3ycn4OHZqfoT3zsQ1b1I5ZUhiVlyfCt/KyiFwAHFTVfck2VXWBB0WkPdvOrfJyfhoaEku20rVHKc4oxZIq6MrLE+FXefkc4GJVfVZEpovIGclqzCIyhcSlpgmQzbUVR2AJp6qvAZ+PrrwsInOAF4G1IqLA70l8J9wvROTnwO0kKjSbPNmK/+gJtIecqfIy8Hdpnv69IGOpNLms+C/Hr4WKKhsDjriJVjG2UchosoSLsGxzZX7JaCv+o8k+7iLM/ywV971k9BuFNMVjZ7gI8ztLZbtktFHIaLKEizC/+9LGU6rORiGjxxIuwvzOUuO5SbS5eZjOzkFOnhyks3PQki0CLOEizO8sZZeMpckSrsiyDftnOkvZJWNpslHKIsr3ywvtJtHSY2e4IrLJ6cpjCVcAE10NYpPTlccSLk/53Dlt5egqjyVcnvK5LLSRxspTMQk30cu+bLJdFibfd+HC+XaLjCmfhPP7xQ6yYI7fZeF43tcmpytLWSRctl/sIEcD/S4LbRTSjFYWCZftFzvI0UC/y0IbhTSjFasQbBPwfRKlFXapame6tvG+T7Zf7PHcqnK6KOrYL4L32waZJ6DtFhkzWuiFYD0PAQ8CDwObfdrGJdvwerbRQL9L0nz6fzYKaUYLvRCsiEwDFqnqp6o6ADSKyIw0beM++2b7xc42Guh3SZpPP8xGIc1oxSgEOwv4OOV5cRLfMTC67WxOV2UeI7VY6YoV8MADNWzZMpMPP5zMvHlDtLYeZ8WKEySftmJF4s/IfST+9iuKmsl4C6Ym37e/v5/a2toR7xsV5VhwNUilVgj2b8DUlOfVAJ+mafP9qUZXmbrzTrjzziGOHPmrt2269ye7bH2tQhVMjXJlrKjGVm5xhV4I1rtkPCgiNSIyFTikqsfTtJ0MMLYR/C5JrR9mCin0QrDe5tVAK/AjYJVPWyj8+lrWDzOFVJRCsKr6DvDOqOeOaQuT371ldt+ZKZSSXfKwYcOGYodgTM4c13WzP8sYUxBlsbTLmFJhCWdMiCzhjAmRJZwxIbKEMyZEJTktENZtPwWM6ybgHhJrRm9WVfXaf0riiyj/pKpXhh2Xt60D+DqJY3NHRI7Xi8CXABdwVbUxXawBxvVNYIOqXjGq/RLgUhInqidU9Ui6Nr99l9wZLszbfgoRl/fd5SdUdSmwFdjotdcDb6lqXcDJlvF4icjXgEe9GJK/wMU+XjOAVaq6kETSPeUTayBU9RVgWppND5D4P3wK7/8xQ1tGJZdwhHjbTyHiUlVXVZ/z2v/I6TsgVgA/EZFfishZAcWUMS7PcuBxEfkvbx1rumMY9vH6RFX/4rV/G/h1ulgDiinViKVF3pk/7v1/9gDL0rVl22kpJlwhbvsJM65U3wK2AXif7IuAl5NtYcelqluARuAjYA3pj2Exj9cy4A8ZYg1baryQOFbp2nyVYsIFettPAHEBICIXAAdVdV+yzftkfBCYElBMWeNS1TiJheONpD+GxTpeMWBIVYcyxBq21HgBBjK0+SrFhIvqbT9p4wIQkXOAi1X1WRGZLiJneH07RGQKiUvNoPjFlbyVfQbwaoZjGPrx8iwHfpt8MDrWgGIaQ0Qmi8gMVe3C+zASkfOB36Vry7a/kku4qN72kykuEZkDvAisFREFfk9iBPAXIvJz4PaU+EOLy9v8qog8DHwXeNxrK+rxSnnKcmB3yuN0sQZCRL4CLBKRL5PoW67zNt0nIq3APwA/9mnLyBYvGxOikjvDGVPKLOGMCZElnDEhsoQzJkSWcMaEyBLOmBBZwhkTIks4Y0JkCVehHMe5yXGcTx3HOcdxnJ2O49QXO6ZKYCtNKpjjOP8BzAXWu67bXex4KoElXAVzHOd84BXgq67r9hU7nkpgl5SV7WoSC5S3FjuQSmEJV6Ecx/kn4DjwAnCN4zg/KHJIFcEuKY0JkZ3hjAmRJZwxIbKEMyZElnDGhMgSzpgQWcIZEyJLOGNCZAlnTIj+H9EWzLD64s1TAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 216x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "train_x = torch.linspace(0., 1., 21)\n",
    "train_y = torch.pow(train_x, 2).mul_(3.7)\n",
    "train_y = train_y.div_(train_y.max())\n",
    "train_y += torch.randn_like(train_y).mul_(0.02)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(3, 2))\n",
    "ax.plot(train_x.numpy(), train_y.numpy(), 'bo')\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('y')\n",
    "ax.legend(['Training data'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The PyroGP model\n",
    "\n",
    "In order to use Pyro with GPyTorch, your model must inherit from `gpytorch.models.PyroGP` (rather than `gpytorch.modelks.ApproximateGP`). The `PyroGP` extends the `ApproximateGP` class and differs in a few key ways:\n",
    "\n",
    "- It adds the `model` and  `guide` functions which are used by Pyro's inference engine.\n",
    "- It's constructor requires two additional arguments beyond the variational strategy:\n",
    "    - `likelihood` - the model's likelihood\n",
    "    - `num_data` - the total amount of training data (required for minibatch SVI training)\n",
    "    - `name_prefix` - a unique identifier for the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class PVGPRegressionModel(gpytorch.models.PyroGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        # Define all the variational stuff\n",
    "        variational_distribution = gpytorch.variational.CholeskyVariationalDistribution(\n",
    "            num_inducing_points=train_y.numel(),\n",
    "        )\n",
    "        variational_strategy = gpytorch.variational.VariationalStrategy(\n",
    "            self, train_x, variational_distribution\n",
    "        )\n",
    "        \n",
    "        # Standard initializtation\n",
    "        super(PVGPRegressionModel, self).__init__(\n",
    "            variational_strategy,\n",
    "            likelihood,\n",
    "            num_data=train_y.numel(),\n",
    "            name_prefix=\"simple_regression_model\"\n",
    "        )\n",
    "        self.likelihood = likelihood\n",
    "        \n",
    "        # Mean, covar\n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.ScaleKernel(\n",
    "            gpytorch.kernels.MaternKernel(nu=1.5)\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean = self.mean_module(x)  # Returns an n_data vec\n",
    "        covar = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean, covar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = PVGPRegressionModel(train_x, train_y, gpytorch.likelihoods.GaussianLikelihood())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performing inference with Pyro\n",
    "\n",
    "Unlike all the other examples in this library, `PyroGP` models use Pyro's inference and optimization classes (rather than the classes provided by PyTorch).\n",
    "\n",
    "If you are unfamiliar with Pyro's inference tools, we recommend checking out the [Pyro SVI tutorial](http://pyro.ai/examples/svi_part_i.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "152d9e4ae9ae421ab0f9f69b8dd67c8b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=200), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "CPU times: user 17.7 s, sys: 460 ms, total: 18.2 s\n",
      "Wall time: 2.75 s\n"
     ]
    }
   ],
   "source": [
    "# this is for running the notebook in our testing framework\n",
    "import os\n",
    "smoke_test = ('CI' in os.environ)\n",
    "num_iter = 2 if smoke_test else 200\n",
    "num_particles = 1 if smoke_test else 256\n",
    "\n",
    "\n",
    "def train(lr=0.01):\n",
    "    optimizer = pyro.optim.Adam({\"lr\": 0.1})\n",
    "    elbo = pyro.infer.Trace_ELBO(num_particles=num_particles, vectorize_particles=True, retain_graph=True)\n",
    "    svi = pyro.infer.SVI(model.model, model.guide, optimizer, elbo)\n",
    "    model.train()\n",
    "\n",
    "    iterator = tqdm.notebook.tqdm(range(num_iter))\n",
    "    for i in iterator:\n",
    "        model.zero_grad()\n",
    "        loss = svi.step(train_x, train_y)\n",
    "        iterator.set_postfix(loss=loss)\n",
    "        \n",
    "%time train()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we are only performing inference over the GP latent function (and its associated hyperparameters). In later examples, we will see that this basic loop also performs inference over any additional latent variables that we define."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making predictions\n",
    "\n",
    "For some problems, we simply want to use Pyro to perform inference over latent variables. However, we can also use the models' (approximate) predictive posterior distribution. Making predictions with a PyroGP model is exactly the same as for standard GPyTorch models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x11e3ffeb8>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ4AAADPCAYAAAAEcj8AAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2deXhb1Zn/P1erbXmR7chx7NhZSdiyQG4IW4CylK1MaUunZFraMp2UDkPptDPTH+0UGHjm90x/z0xnpkxXKDOkGUpaoOyENaWsTbgJSQhJyGo7jjd5kSVrvVf3/P64dmInkm3Fki3b5/M8AUv36uiVLX117jnv+34VIQQSiUSSCbaJDkAikUw+pHBIJJKMkcIhkUgyRgqHRCLJGCkcEokkY6RwSCSSjHFMdAAjcdddd8n9YolkgvjhD3+opLo/74UD4L777hvxnI6ODqqqqsYhmlMn32PM9/gg/2PM9/hg9DHee++9aY/JSxWJRJIxUjgkEknGTIpLlcEYhkFbWxvxeJzB6fKmaRIMBicwspEZS4yKouDz+SgtLc1yVBJJ5kw64Whra8Pj8VBbW4uiHF+30XUdp9M5gZGNzFhijMVitLe3S+GQ5AWT7lIlHo/j9XqHiMZ0wO12YxjGRIchkQCTUDiEEKMSjQ0bbCxa5KKw0MWiRS42bJh0L3UIiqIgK5klY+VoIMqm/YExjzO5P01p2LDBxu23O2hqUhBCoalJ4fbbHRmLxw9+8APWr1/PsmXLWL9+Pd/97neHPX+47avB3H///bz55psZxSKRjJWjgSi/ePMwe9ojYx5rSgrHPfc4iESGzkoiEYV77slsSec73/kOt9xyC1VVVdxyyy1873vfG/b80eSbAMydO/ek2YNpmvz0pz/NKD6JZLS0BGI89HYjANmYt066xdHRcORIZveno6KiYsjt8vJy7rvvPoQQFBUVEQqFOPPMM+nt7eX666/nBz/4AT/84Q+5/fbbue6663j11Vf53e9+d+zxL7/8Mrqu8+qrr/JXf/VX/PznP6e6upr33nuPO+64g+eff57Vq1fzwgsvcNZZZ7F9+3buueeeTF++RDKEtt4YD77dgF2B4gIHwb7EmMfMqXCoqnoJcK+maVeccP+FwEVYM57/0TStI5vPW1cHTU2p7x8r8+bNo7a2liuuuIKtW7cSCATYtGkT3/jGN3A6ncyaNYuKigrWrl3L448/PuSxjz76KI899tixLdlLL72UQCDAoUOHqK+vp7a2lqVLlwLQ09PD/v37xx6wZFrT1hvjl29ZouEtchLXk1kZN6eXKpqmvQkUpjj0L8C/AY8Bo5vfZ8D99xsUFQ2dkBUVCe6/Pzu7Eg6HAyEEv/71r7nssstOuuwYWLw98f5Dhw4NWdx97LHHWLRoEV6vd8h569atY+XKlbjd7qzEK5metAetmYbSLxrZZDzWOIbMi1RVXQQYmqYJTdOagNXZfsKbbzb52c8M6usFiiKorxf87GcGN99sZjzW/v37OXz4MAcPHgTgwIED7Ny5E4Ddu3ezfv16/H4/fr+fAwcO0NjYyKFDh2hqaqK7u/vY4wC+8IUv8M1vfpOdO3dy8OBBAoEADz/8MN3d3ezbtw/TNNm8eTPNzc088sgjdHZ2cvTo0ez8UiTTivZgnAffagAE5VkWDQAl11t8qqq+oWnaZYNuXwh8R9O0m/pvH9U0rTbd4++66y7xN3/zN8duBwIB5s2bd9J5yWQSu92ezdCzzlhjPHz48Emzk2wSCARyOn42yPcY8yG+zrDO+q0dCCHwFg5djYgbJpFonL+74uTP0In89Kc/zavqWD9QNOh2fKQHDK7kCwaDabMv8z1zFMYWo81my3nlZb5XdkL+xziR8flDcX6vHaawoJByz9D32rZNpbz4SBW9fidP1VuX9KcyC4dxFA5VVe1AkaZp+1VVLei/bz7wxnjFIJFMZTr74jz4dgOmSUrReOKBGvS4tTrR1AS33+4ATk08crrGoarqEmCBqqpnA9cBP+g/9M+qqn4X+BLwj7mMQSKZDnT1JfjlWw0YSXGSaABsXFd1TDQGOJXcpgFyOuPQNO1DYGATdBfwXP/9m4BNuXxuiWS60B1O8ODbhzEMQUVx6kvhgN+6v3hpE3qPh/iRSiDz3KYBpmQCmEQyXegOJ3jwrQbiukllsSvteV5fArGwCe/qfZhxB0cfvAwz4j7l3KYpmXKeLR5//HEuvfRSHnjggRHrVOB4DcqOHTt46qmn0p7X1tbGz3/+82yGKpmG9EQSPPR2A9FEcljRMIVg0Vc/wLt6H8KEnjdOx4y4x5TbNGlnHIvvfS1rY31835Up71+5ciWnnXYad955JzfeeCM7duxg2bJlaccZqEFZtmxZyvN+8pOf8O1vf5vq6mr++q//OmvxS6YfgYjOQ281EkkkmVGSXjSMpMkre/y0KxEUoRD+w1L6dtRSXy8mx67KZCcWi/EP//APLFmyhM997nPs27cPTdP41re+xcGDB4lGo8dqUB5++GEKCgpYvXo17777Lu+88w5r167l5ZdfZvny5SQSCbZv384dd9zB+vXrmTFjBkePHuXqq6/mvvvuY8WKFTQ2NvJf//VfE/2yJXlIb1Tnobcb6EsY+IaZacT0JC982E5rMI7bYeP6s2dSeXEPwb6j3P/Z5WOKYdIKx4mzhFx1AOvo6OCJJ57gtttu45lnnuFHP/oR3//+97n00ktRVZXDhw+zfv16Hn300WM1KAsXLqSxsZFf/OIX3H333VxzzTUoikJNTQ2f+MQnaGhoYMuWLTz55JPMmTOHa6+9lk996lPceuuteDwe/v7v/55PfvKTWX8tkslPMKrzq7cbCcYMqoaZaYRiBs/tbKM7olPstnPD0moqPa6s1apMWuEYL6qqqrjpppsAeOGFFwArA9Tr9XL11VcfS0dP1WDINE127959bAZx4vFYLIbf7wegpqaGwsLjZT2yaY/kRIJRnYfeaSQQTVBVkr6OqSuc4NkdbYQTSSqKnNywtJqSgux+1KVwDMOWLVs4fPgwXV1d6LrOxx9/TGdnJ2vXruXLX/4yZ5xxBvfffz9r1qzhm9/8JgUFBRiGgRCCpqYm1q5dy9e+9jXOOecc7r77bhYsWMCvf/1rampqaGho4M477+T73/8+Ho+H66+/noaGBlpbW2lsbKSlpQW/34/P55voX4MkDwjFDB5+p5GecIKZpelFoyUQ44Vd7cQNk1llbq4/eyYFzuyXYuS8VmWs3HXXXWJwg5wDBw6wcOHCk86b6s2KIf1rzxZTyUxooshFfH0xg1+900hXOJ5yprFtUykb11UR93bhu2E7isNk/owiPnmGD4d96MZpXE8S7IuMao3j3nvvTVurIrdjJZI8pi9u8N/vNtLZl140nnigBqOmFd+nt6E4TMI765jZsfAk0cgmUjgkkjwlmkjyyHtNtAfjaS9PNq6ronDZYSqv3oVig8Bbi+jcuISX1s3MaWyTbo1joNv3dLNHmI6veToT15Os39xESyDKrLKCtOeZ85oov+xjALpeOpu+HXOA4ynmuWLSzTjcbjeBQGDa7TrE43Ecjkmn85JTIGGY/Ob9Zg53RqgeZiF0a2OA8ss+RgjofHHpMdEA8Pr0nMY46d6J1dXVtLW10d3dfZIFpM2W3zo4lhgHLCAlUxsjafLbrc3sbQtR6y1IO8vUGgP86XAPCAi8uoTwh8eLTpxuk2u/ktU2vicx6YTD4XAwe/bsk+7P99V2mBwxSiaOpCl48oMWPmwOMrt8FKIBXHHGDKLCxsajCQJ+J16fzrVf6eDcy1N7FAdjSQqdY/+CnXTCIZFMRUxT8OyOVrTGwLCi8X5DD5sbLCe2K0+fwenVJVAdTCsUx8YXgrZea2fm2oVj9x+WwiGRTDBCCDZ+1M67h7qZ7S3AlkY0tjT0sKUhgAJcebqPxdXFoxo/YZi0BeOsqPdy4/JZBHu6xhyzFA6JZAIRQvDaXj9v7Ouk1luAzZZaNDYf7uH9xn7ROMPH4pmjE41QzKA3pnPj8llcOL8CRVEYfm4yOqRwSCQThBCCN/d38cruDmrK3NhTiIYQgi0NgYxFQwiBvy+B027jG6vnMW+GJ6uxS+GQSCaIPx3u5vkP25hV5k6Z5bl1Uwl/2BmiYHkAYcIi6lk8c+S6k6QpaAvGqa8o5C9W1mXdjAlybwH5d0AHUKZp2k8G3f8ZoLL/ZkTTtN/kMg6JJN/Y1hTg9x+0Ul3qxplCNLZtKuWVLWFKzj+MMBU6n1tO6+FqZtAy7EJoXE/SHkpw8cJKrj1rJi5HblIUcpb4oKrqxUClpmnrgXJVVVcNOvwtTdN+pWnar4C/zFUMEkk+8uHRXn6rNTOzxJX2g/36tjAl5x9AmND53HIiey1rg43r0m/n90Z1uiMGa1bO5s+WVudMNCC3M47rgD39P+/uv725//ZWVVXvB54FfjbSQB0dIyezBAKBU4tyHMn3GPM9Psj/GEeK70BnlN9u76S8yI6RMDFSGMfvbo9QqPYhBHS9sIzI3prj4/udhMPhkx7TFTYocNr4wtIZzCpIHOvzcioxjoZcCscMoKf/5xhQPejY3cBDwL8Cnx9poNEmTU2G5Kp8jzHf44P8jzFdfAf9YZ7b56emsoQiV+q1io9agmxu6gOg++UlhHcPTXb0+nQ8nqELne3BONUVRXztojmjXs8Y6+8wlznag60eS4DBm8f/F/gGlmP9hhzGIJHkBY1dEf773UZKCxxpRePjtj7+sM/6mMwza4jvHSoaqVLJ24NxKj0u1l48etHIBrkUjheBpf0/nwm8rKpqWf/tpZqmhTRNewHI7+47EskYORqI8qt3GvC47BS7U0/yD3SEeW2vdXlxwfxyrr/czU13tuCtSoAi8FYluOnOoQujA5mgay+eS2nh+H6McnapomnaO6qqfkJV1VuBQP+/XwBrgP9QVfWbwFHgl7mKQSKZaNqDMR56uxG3w5a27+fhzgiv7OlAACvneFlRb7ndn3t56lRyIazt1hpvIbdeUI8njRjlklxbQP7zCXet6b//+Vw+r0SSD3T2xXnwrQZsCpSlmRE0dUfZ+FE7poBz6so4b6532DGFELT2xplTUcSXL6hPe9mTa2QCmESSA7rDlsuaKaAihQk0WJcwL+6yRGNJbSkXzi8ftlnTgGjM93m4ZVVdTpoQj5b8bmAhkUxCevu9T2IJM61otPXGeP7DdgxTcGZ1MZcsrBhRNFp64yysKubLEywaIGccEklWCSeS/O87jYTietrmwi8/V4Dn6q3YCgQzhJfLFntHFI2jgRhnziphzcq6nCZ2jZaJj0AimSL0xQ02bPfTHU5tmLRtUylPP1pC0VXbsBUYhD+uZucDq9j+h7IUo1kMiMbZtaV5IxoghUMiyQp9/YZJnWE9bUfyl5/2UPnZLdiLdCIHquh89hz0mCNtGrkQguaeGMvqyrhZnZ03ogHyUkUiGTOhmMGv3mmgqy+BL82aRm9Up/Cqrdg9CaKHZ+B/+lwwLSFI1ZHcmmnEWTq7jD8/tzanHimnQn5FI5FMMoL9zvFdfemtGYMxnae3t+IoiRNrqsD/exWSxxc3T+xIPiAaZ9WU8AU1/0QDpHBIJKdMb1TnobeH93Ptixs8s72NUDxJiSii57kVCOO4aJyYRj6we3J6dTE3q7NTltznA/JSRSI5BQa2XAPRBFVpRCOSSPL09jZ6Ywa+Yhc3LpvBR6KDjeuqUnYkH8jTWOjz8Bcr82tN40SkcEgkGRKIWJcnwZiRcvcELPvGp3e0EojqVHpc/NmyatxO+4hp5PNmFPGlVXW4JzhPYySkcEgkGdATSfDQW430JQyqSlwpz4nrSZ7d2UZ3WKe8yMmnl1VTOIIQtAcT1JUXccuq+glP7hoNUjgkklFipZE3EklYlx6p0JMmL+5sx9+XoKzQwaeXVY9YT9IejDPLW8BXLqincIJqTzJFCodEMgq6wwkefKuBaCLJjBQzjW2bStn4aCXuyz6goC6OWzi5cVl12jL6AdpDVmn8VyewYO1UyN/VF4kkT+gIxfnFm4eJ6elF48mfVeG+ZAcFdT0YoQKa/udC9r9bMfy4wTiVRS7+8sI5IwpMviGFQyIZhuaeKD//42GMpKAyzeXJxkcrKP+zrRTM6cIIuWnfsIqYv3jYxsL+vgTeIhdfu2gOxWn6dOQzky9iiWScONQZ5n/eaaTAaUvbYSthmLgv294/03DT/tj5GD2WYVKqjFCAzr4ExW4HX7tozrh37soWUjgkkhTsbgmyfssRSgscaS8jEobJczvbKKiLY4QK+kXjeCPhEzNCwRKNIpdj3HuEZhspHBLJCWxtDPDbrUeZ4XGm3eVIGCbP7myjLRjHJZy0P3HeENFI1Vi4q0+n0Gln7cVzKC9KfdkzWZgQJ7f+Y6cDq4Fdmqa9l8s4JJLRIITgrQNdPP9hG1UlbtxpMjfj/aLRHoxT4rZz4/KZHBS9bFznTpkRCtAd1nE5FNZePJcKz+QWDcihcAxycvuRqqp3q6q6StO0zf3HFgNrNU37u1w9v0SSCaYpeHVvB6/v8VNdltqWEQaSu9ppD8UpcTv4zPJqSgudxzJCw+HwSb4nPWEdu80SjXQLrJONXO6qpHJyG+ABoFFV1R/3C4xEMmEkTcEzO1p5fY+fGm9BWtGI6Ume2dlGeyhOaYGDz5xTPeLiZk8kgaLA11fPxZcmPX0yMu5ObqqqeoC5wE+A2cBmVVXnaJqWwgzPQlpAjg/5Hh9kP0Y9afL87m4+ao9QXeIkFo2kPC9umLz8cYCuiEGJ28bVi8qwJxOEw0PftrFY7NjPvTEDBHxpxUyI9tIRzWrop0y+W0Cmc3JzAVFN00ygSVXVFixRaUo3kLSAHD/yPT7IXozRRJLfvN/MoV6ThdXp+35G9SSv7GmjK2JQVujgM8tmDZt74fF4CER03AUuvrF6LtVlBVmJN5vkswVkSic3TdN6gLiqqsX9x/xYxkwSybgRiOg8+HYDB/191HoL0opGOG7w9PZWOvsSeAudfGb58KIBVsm9YQq+fnF+ikY2yJlwaJr2DhBL4eQGcAdwr6qqa4D/p2laMldxSCQn0tob46d/PERnX4JZZelFIxjV+f0HrXT1V7neuHzk2pNQPEnCMFl78VxqvFNTNGDinNzeB97P5XNLJKk46A/zyHuNuOy2tGXxAF3hBM/saCOSSGJ0lrDzsVU0FCsnbbMOJhjVSRiCOy+fx+zywly9hLxAJoBJpg3bmgL8VjtKeZFjWL/VtmCM53a2EzdM4s0VtD+uIhJOAhF44oEagJPEIxjViRsmX1xRRd0UFw2QRW6SaYBpCl7f28EGrRlfsXNY0TjSE+WZ7W3EDRO9aQbtvz0PkTi+5arHbScVr/VGdeJJwddXz6OmdGrkaYyEFA7JlEZPmjy1vZWXPupgVql72JZ8B/1hntvZhm4KFs/00PI7dUhj4QEGF6/1RnX0pOC21XOn/OXJYKRwSKYskUSSX//pCFsaepjtLRjWZmBPa4iXPurAFLC0tpQrT/fhrUy9Zj9QvBaIWLsnt62eR613+ogGyDUOyRSlJ5Jg3XtH6AjFqPW6T9o52bap9Fi38apLD1C4qhOAlXO8nDfXyum49isdPPFADXr8uOAMFK8FIjpJAV+/eN6U3j1JhxQOyZTjaCDKI+81ETdMZqXIo9i2qbRfEBTKLt5H4aoDAMw1a1g173ha+MAC6Il2BvPP78IUcNvquSnHnw5I4ZBMKXY297JBO0qRy5a2ofDGdVXoCYXyKz+idEUjwoSujUvp7ZjJpy4/MOTcE+0MeiI6Aks0pmpy12iQwiGZEiT7d05e2+vHV+wa1mKgt8eG78atFC1qRxg2/M+eQ3R/NShi2OfoiSQAhdtWz2Vm6fQVDRilcCiK8q/AvwghunMcj0SSMZFEkie2HWXX0SA13gLsttSZoAPn1tzyHg5fkGTMgf/3KvEjlUDqjl0D9IQTKIrC11fPS2v3OJ0Y7YzjMeBqRVF8wF7gVSHE8PIskYwDnX1x1r13hK5wnNnl6dPHwZoxPLezHYfPwAgW0v67lRhdJUDqjl0DdPXp2O0KX79YisYAoxWOD4FdwK1YvTSeVBTlCPCoECKUq+AkkuE46A+z7k9N2FBGXKRsCcR4YZeVDeordjHHnMfrdjcBRaTs2DVARyhBsdvOX100dZrwZIPRCsdrQCHWzEMVQvQpilIIvABcnqvgJJJUCCF491A3z+5ow1uUvpnwAPs7+nhtTydJIZhbWcjVZ1bhtEc5/4SF0BNpC8ap9Fi+J5O5sXAuGK1wvAX8kxDCGHSfAezIfkgSSXr0pMnvP2jlT4e7qS51D+voLoTggyO9vHvI6id1dk0JlyysxDbMGsjA41qDcWq9hXz1/PpJ6XuSa0b1GxFC/CDFfTrw7axHJJGkoTeq89gHfjpjNmZ7C4YVANO0Gg9/2GJdSV84v5xz6sqGXQMBSzSOBuIsrPLwpfPqJo2X63gjpVQyKTjoD/PoliOE+nTqfKXDCoCeNHl5dwcNXVFsClx1ho/TqorTnj+AKQQtgRhnziplzcrZw85mpjtSOCR5TdIUvLHPzyu7O/AWOan0OIYVjfc2FbK5qxmHL4oZd3C2cw6nVZkjPo9pCpoDMVbUe/ncuTVpGxZLLKRwSPKWYFTn8a0tfNweYlaZG4fdRtiIpz3/j5vsbI8cxOGLowcK6Xj8PNrCRZTSkrb5Dlji1BKIcdGCCm5YOmvYPBCJhRQOSV5yuDPM/24+QjxpjpifAbC3LcROsxt7sUmsqQL/0yswo9b26cZ1VWmFQ0+atPbGuWKxj0+eWTXiwqnEQgqHJK9I9i9qbtzVTlmhg5lFwydcmaa1Nbu9OYjigNC2OXS/fiaYxy810pk/JwyTtmCc65dUc+lplSOKk+Q4E2YB2X/8YWC9pmlv5DIOyeSgL2bwxAdH2d0SGtZNbYCYnuTl3X6O9FiLoOF3zqD7rfknnZcqlTwcN+iJ6Hz+3BpWzi2XopEhOVsBGmQBuR4oV1V11QnHbwBGXuqWTAuauiP8+A8H2d/Rx+zy9G5qA3SHEzy+rYUjPVEKnDY+vayaK8534XQPXQhNlUoeiOj0xZN87aI5nDevQorGKZDLGUcqC8gB79h5/c+9J/VDhyKd3MaHiYgvaQrePxJi04FePC4bJW47kcjJbmo736zg9d/MprfTReXyFko/uRNTMakodHDFaWUUO028q1q54ba4dV6Xi7LKBFf8RTOLV3UTDlvjdEUMXHaFm5f78CpROrJsr5bvf2PIfye3dBaQDuBaTdN+pqrqstEMJJ3cxo/xjK+rL8Hvtx3lkD9G3YzStHkT2zaV8twvrcY7pasO4bl0L6YClaKMm1TvkNnJBdfGueDagyeM4EEIQVswTk1FEV+9sJ7yotzVneT73xjGHuNEWEBeAnxJVdU/x/KQvVFV1es1TZNubtOEgVTwp7a3oCjKiLsmG9dVYSQFMz61A89ZLQAE3lpE74F6nJ84USROxkrsirOoqpg1582mSGaDjplcCseLwLXA7xhqAbkJuBBAVdV/At6QojF96IsZPL2jlR3NvVSVDN9w59hjknGqv/wBLl8IM26n84Xlo2q8A2D0b7eeN7ecTy+fJRO7ssREWUBKpiEft4X4j9cPsKc1RF15wahEY39HHzVfeQeXL4Te5aHtfy+yRIPhG+8AxPUkLb1xrj6ris+eI7NBs8mEWEAOOv5PuXx+SX5gbZt28M7BLsqLnFSPkJsB1qLp2/1FaooTIh/PovPFpYiE9ZYdrvEOQF/cIBAxWLNyNufWe7P2WiQWMgFMklOO9ETZ8H4zXeEENWXDt/UbIBjVeWl3Bx2hBDYFVi+sJGH6eMlrEvAP33gHrDZ/uglfXz2XBT5Ptl+SBCkckhyRMEze3N/Ja3v9eFx2akfpPXK4M8Jre/3EDZOSAgfXnFllteurDbHi8hDhcBiPJ7UYCCFoDyYoK3Ry2wV1076hcC6RwiHJOoc7wzyxrYWucIKqEtewawvHjJE67VRfsxf3knYA5lUWccXpM0a1DgLWpU1rb4xFM0v4glo7YlcwydiQv11J1gjHDV7Z3cF7h7opLXSMOMsYMEYynQlm3rwVd103wlSYSzXXnX2y+1o64nqS9lCCS0+r5JqzZg5r9SjJDlI4JGNGCMFHLUGe2t5KRE+OaFEwwMZ1VdhrOpl5/Q7sngRGyE3ns+cSTBRzwwj9QAcIxnRCsSRrVs4eVYcvSXaQwiEZEz2RBM/ubOOjliAVRU5qika3rqAnTZRl+5h5biMA0YZKOp87BzPiJj6K/AywOpC77Aq3XzqP+oqikR8gyRpSOCSnRNIUvN/QzQu72hECZntH7pkxQEcozqt7/JSeqyOSCoE3FxN8fz4I6/Ej5WeYQtDaG2d2eSFfPG92TtPHJamRwiHJmNbeGE9+0EJTV4SqEhfuUS5gmkKwramXLQ09mAIKhZvGDSuINpcfO2ek/Aw9KejqiaHO8fKZ5TWyL+gEIYVDMmr6Ygav7/Xz7qFuCp22UXXmGiAY1Xl1r5/WXqv139LaUi6cX85OEWbjOs8QN/h0+RnhuIE/rHPTylouXjCyzYEkd0jhkIyInjTZ0tDDyx91YJgms8rco+7LKYRgb3sfb+7vQk8Kilx2rjx9xrE1iRPd4NON0RnWUYA1y31ccNqMsb4kyRiRwiFJixCCfR1hntnRSldfHF+xC7dz9OsJUT3JG/s6Oei3+mssmFHEZYtnUDjKSxuw1lLagnHqygtZs3I2Rjj/+11MB6RwSFLSHozx/IftfNweoqzQwezywlE/dtumUl59XcF9/h4cJXHswsZlZ1Rw+szijLZLo3oSfyjB6oVWfobLYaMjfCqvRpJtpHBIhtAXM9i0z8+7B7tx220Z7ZYA/GlTIW/u78BzVSsAseZyel9ZSvSrAZTq4S9JBtMd1tGTJresqmNJ7fAGTJLxRwqHBLBqS7TmEJuPdqEnTapLR7+OAdZlzZ62PrboRyhabGAm7ATeXkRImwdCYeM614hrGWDtvLT1xvGVuPnieXVWnYok75DCMc3RkyY7m4O8tLudjp4Qs2eUUpDBOgZYzX/f2NdJcyCGzQGE4d4AABIlSURBVA3RQz66Xzkbo/d4UlY6i4LBxA2T9mAcdY6XTy+bNeo6Fcn4I4VjmmIkTXa3hnhxVzuBqE55kZNZpaPryDVA0hRsP9LLlsYASVNQ4LTR/dqZ+P9UDwydrYyU1BWI6ET0JDedW8N50q4g75HCMc0wTcGethAbP2qnM5SgrOh4MVp4+M/2ENqDcf7wcSed4QQAi2cWc/GCCvbE7TzxgUCPH//gD5fUZQpBezBOaYGTv7lwTkaLsJKJQwrHNME0Bfs6+nhxVzvtwThlhQ5qyzPvV5EwTDY39LCzOYgASgscXLaockheBlgFbCMldQ3smqyo93LD0mo8shR+0jAhTm6qqq4B/hYoBW7RNE3LZRzTGSEEB/1hNu5qpzkQo7TAzuxTEAwhBB+39/HeoR7CiSQKsHx2KavmlZ/Ub2OkpK6BhC4hBF88bzbLZsuq1snGuDu5qaqqABFN01YB/wbcl6sYpjPWomcv//n6QR56u4FAVKfW66a0cORFyhNp7Y3xxLZWXtvbSTiRxFfs4vMrarh4YWXGDYD1pElzT4yasgK+fcVCltd5pWhMQsbdyU3TNAE803//+8AFOYxh2hGOG3xwJMAf9nUSjhmUFjqpzTAXY4BQzODdQ93s78+6cgoHvW8vRnuvngM+Y9i6klQEIjrhRJLrllSzemFlRtu9kvxi3J3cTuBK4N9HGkhaQI5MT8Rga3OIrUfDJE0Tb6EDr8sGyTgpHBVTEovFAKsCdVdbmA/bIiRNsCvgS1TywS/PIRG28ioCHS6e+PEs4vE4Sy/pHnZc0xR0hA3KCx2sWVLJrFKTrk7/Kb3OfP8753t8kP8WkOmc3ABQVXUh0Khp2u6RBpIWkKkRQnCkJ8qb+zvZ1RLCpkBtZcmxy4dj/TxHUXk6MF5zSPDuoW7CiSQAC30eLpxfzn994wwS4aH5HXrCzqYNdVxwbTztmOG4QSCic9kZPq45a2ZWcjPy/e+c7/FBfltApnNy61VVdSawTNO0J1VVLQaEpmmyCmGUxPQkH7eFeOtAN0d6orgdCrNK3UPKzAf6eepxS0QCHS6eeKAGIKV4tPbGeHNfD/6wAYCv2MXqhZXU9G/VpkvgSne/KUR/hy4bf3lBPYurS+RaxhQiZ8Khado7qqp+4kQnN1VV7wBeBgxVVb+HlSmk5iqOqYJpCpp6omxtCvBBUwAjaVJc4KDWm7qp78Z1VcdEYwA9bmPjuqohwtEejPPq1hABJWQ9T8TFooKZXL3CPmRcr08n0HFyRmmqxK6+uEF3WGd5XRk3LKk+pQVZSX4zUU5uy3P5vFOJ7nCCD48GefdQN71RHafNRqXHOWIn75FmCP5QnM0NPTR0RUEBM24nqM0juHkBbTYbVbQMEZhrv9IxZAYDJyd2JU0rmcvjdvC1C+UsYyojM27ykLieZL8/zHuHujnoD6MA5UXOUZsaQfoZQsWibl7c1c6hTmvFVOg2glvnEdwyHzNqna/DSTOTkRK7eqM6oZjBRQsqufKMKukIP8WRwpEn9MUMGrojfNQSZFdLCD1pUuSyU1M2en+RwZw4Q3BUhqi4ZB+Fi9o41Al2m8KSmhKeu2slZvhkQUo1Y0mV2KUnTTqCCWaUuPjSqnrmVspu49MBKRwThGkK2kNxDvr72N4c5GhPFAG4HTbKixxjdlYf+IC//JQH5fQGPGe2oChgU+DsmlJW1JfhcTv4o8dGIMWy9EhFaUIIusMGMSPJVWf6WL1whmwcPI2QwjGOxA2Tfe197GkLsbM5SDhhoADFBXZmneLMIhWi3z6g1ddOyeesSxKbAmfOKkGt91JccPzPPpq1i5Neh56kI5RgbmURnz2nhuoy6dE63ZDCkSOSpqArnMAfinOkJ8qhzjAHWntwuy2XM2+RE+8ozYtGi2kKDvjDbG/upSNkVa3aFDi9uhi13ptydyOTorSkaW2x2m3w2XNmsXJuhcz+nKZMGeE41BUjbAtTWuCkpMAxrtPmpCnoDifoCMVp7olyqDPC0UAUUwgQoNjA43Iwo8hBSUn2v50ThslHrSF2NvcSiluJWwUOG0tqS1lSWzriQuXA2kU6J3ghBF1hnZiR5Px5FVyx2Ce3WKc5U0I4jKTJb7f7KSzsAwSmgGK3gxnFLnzFLqpK3JR7XJQWOCh2O3DZbdhsCnYFFEXBblOw9f88GCEEMd0kqicJJ5JEEkkiCYNg1KAnotMbtf51hOKYQiAE2AZEoth10rdxOJmdb+eBjNBQzMC3+iBFZzWTVEwAvIVOlteVsnhm8ZjXScCqVwlEdU7zFfOppdXMkpclEqaIcAAIBLPKrDoKIQR60poFtPbGSBgmKEN7Ugkx8DiOPcZht2FTwG6zYVcgZpiYQqCgoCggTIGJNY7TruC023DalZQikSu2birh+aedFJy/i9rFbSg2QRIoFR5WLylmbmVhVtZK4oZJZyhOucctMz8lJzFlhGMwiqLgcii4HDZKRvkYIQQCS1AGfi4tcOSNW1ggorOnLcT70RYqP2vVhghTIfxRDcH359MrCpm3bnQO78NhrWPEcdhs3LBsFivnlMvdEslJTEnhOBUURbFmJMf/M+EkDJMD/jB7WkO0Bi2xsHlA7ykivGs2fR/OJhmyWu0lRunwng4hBN0RA5se58L5FXxisY+SAvn2kKRGvjPGgUyqVIUQHA3E2NPWx0F/GMO0BMFpU1hQ5UH778V0fTiTTJsBp+P4wqdJndfN51ctkNurkhGRwpFjRlOlaiRNjgZiNHRHaegMH9sZAagpK+CM6mIW+Dy4HDYqrjF4Yt/omwGn45hg6CaLZxZz1Rk+XHqImVI0JKNACkeOSVel+tITpbgWN9PQHaG5J3ZsZgFQ4rZzenUJp1cXU3bCtmcmeRepsDI+daJGkkVVxVx1RhV15daCakdH3xhfrWS6MG2EI9OmNtniWM2HYuKuCVC4oIPCBR24qkK8sf/4eb5iF3Mqi5hTUUh16fBZpKNxeD8Raw1DJ5I4Lhj1FdnZgZFMP6aFcGTa1GasCCEIxQw6+hLMvHo3orQP16wA9gLj+Dm6nQWz3MypKGJOZSHFbgfbNpXy3/+UXXETQtDTLxgLfMV88swq5kjBkIyRSS8cGzbYuPtuJ0eOXIi3KvWHbbRNbU6FYyIRStDRF8cfiuMPJYgZVkJWwbJBz9ldRPRQFXqTj+tvjKGeffzSINviljQFnX0J9KRggc/DVWf4mFtZJAVDkhUmtXBs2GDj9tsdRBNWWla6D1umbe8GGMgcDSeShOPGsf/3xZPHfg7GDOL9IjGYAqeNqmI3vhIXfQ1e3n+ijp6GUrw+nRu+0sG5lw9dT8iWuCUMk84+q07lnLoyLl44I6sFdBIJTHLhuOceB5GYYM4/vASAqdsQup23ozZ2bzZx2BUcNoWaL/aiRxwIw44wFRSbQLGZOIuSPLMjjCmE9c+0emUm+zNPw3EDcxTpEYVOG74SN1XFLuv/JW6K3YNa782Hqy4/Sji8L2UtCJy6uA0QSSTpiei4HTYuP93HyjnleItkPYkkN0xq4ThyBBSnianbsDlNbE4TnNa3fyB6/DxnbZx0H6EjPWkO9GMXNko9dorddjwuBx63HY/bQfvHJWx9oYrAkWJKCm2c8xU/5y499cueTHp6DiCEoDdqzYC8RU4+s3wWy2aXSZd3Sc6ZKAvIRcAXgAjwnKZp+05l/Lo6aGpycOTfrwUEijOJ4khSVh3jr//1MLopMJICwzT5eHshH7xVTKTPRpFHcM4lQRYujWJTFGz9RW4HthWzaUMVesyOMOwk+9w4bDZuurOFc887LgrbNpXyh8HrESHGvNiaSV8MPWnSHdYxTEFdeSGfX+HjtKpiWeIuGTdyJhyDLCB/pKrq3aqqrtI0bXP/4R8Dn8dqb/kY8NlTeY777ze4/XYHkYgCKAjdgcNm45qb/JR7hn57118OV10+uNWVDRh62fDIunlETvjWT9V/MxeLrSPlZ5hC0BvViSRMnHYbK+eWc269l9mn6NImkYyFcbeAVFW1EFigaVofgKqq81RVdWiaZqQZJy0332wCBnffbefIESXtrspoGe06w1jXI9KRKj8jkkgSiOgILHOkC+ZXcFpVsSw8k0woE2EBWQ4M/nQYgA9oTTfQcBaQl18Ol1wm+OeXDlE/w5pBhE/R2qmsMkFvpzvl/eFBg472vBMZsFgcCcMUBCJJDGFSVuDgwtklnD6zkLICBxAj0D26cTJlutgX5pJ8jw8mrwVkFzC4IKIIy6wpLSPZ1RlJE5f7SNodi9Fy3a3+lOsM193qHzL2aM9LRbrjMT1Jb9TAMAVOu42LF5dx7hwvs72F41raPx3sC3NNvscHk9cCslFV1SLABI5omhYdbqDxYrR1IGOtFwFrzaIvZu2IAHjcDlbNK+f06hLqywtxy50RSR4z7haQWG5u/wf4LhAHvpOrGE6F0daBnEq9iJ408Yfi6EmrUdCcikI+sbiMBT4PVSUySUsyeZgQC0hN03YBu3L53BNN0hRE+rNLTQGKAkndRJ1fxlk1JcypKMLjntRpNJJpjHznZgEjaVoikUhavUwVsCsKtd5Czq33UldeiK/ETTIcYFb1zIkOVyIZM1I4RkAIgWEK4oaJbpjEkyaGKbAP6n5st9mo8xZwns9DrbfA6qpedHID446ovBSRTA2mhXCIfuuCpBAkzeP/jEH/N03rcsLSg4EPuGW1UOS0U+5xUlnhotLjorLYqkUp6bdbKHHnT1NjiWQ8mDLCoSgKrb0xjvkg9Psf9HsiYbcpFDhsuBw2itwOipx2q+5koP7E5aDAacfdf47bYaPIZYlDNvxJJJKpxJQQDofdxprlPsorKnDZbTjtNlwOy/fE1e994pAffokka0wJ4QCYW1FAVVXxRIchkUwL5NewRCLJGCkcEokkY6RwSCSSjJHCIZFIMkYKh0QiyZhJsaty7733TnQIEolkEIoQY3M5l0gk0w95qSKRSDJGCodEIskYKRwSiSRjpHBIJJKMkcIhkUgyZlJsx55Irh3issEwMa4B/hYoBW7RNE3Lp/gGHX8YWK9p2hvjHdugGNLGqKrq6cBqYJemae/lU3yqqn4GqOy/GdE07TcTEV9/LJcA92qadsUJ918IXIQ1efgfTdPSe5CkYNLNOAY5xK0HylVVXTXo8I+B/wB+AvxwIuKD9DGqqqpgvZFWAf8G3JdP8Q06fgMwoaXGw8WoqupiYK2maQ9NoGgM9zv8lqZpv9I07VfAX05EfANomvYmUJji0L9gvQcf4xTeh5NOOEjtEMdghzhN0+LAPFVVJ2pGlTJGTdOEpmnP9N//PsOYUOWYlPGB5ayHNRPdk+Jx40naGIEHgEZVVX/c/wGeCIaLb6uqqverqqoCPxv3yE4mMfhG/8zc6H8/NmHN3DJiMgpHpg5xE0G6GAdzJfDv4xbRUFLG1y+012qa9tQExTWYdDF6gLlYs8ofAY+rqupKNcBExNfP3cAC4F+BN8c5rtEwOHawPjsZMRmFI2sOcTkkXYwAqKq6EGjUNG33eAfWT7r4LgG+pKrqG8BXgf9UVbV23KOzSBejC4hqmmb2f1u2kFqYJyo+gP8LfAPrUmDDOMc1GgbHDpa/UUZMRuF4EVja//Ngh7g41vS1SFXVAibWIS5ljACqqs4Elmma9qSqqsX936B5EZ+maZs0TbtQ07TLgEeAv9U07egExDdcjD1AXFXVgTUYPzARMab9GwNLNU0LaZr2AjA2J/IsoqqqXVXVEk3T9tP/Jauq6nzgjUzHmnTCoWnaO0AshUMcHHeI+zYT6BCXLkZVVSuBl4HvqaqqAX/E2gHKi/jGO47hGCHGO4B7+3eo/p+mack8i+8/VFX9pqqqnwV+Od6xDUZV1SXAAlVVz8Zah/lB/6F/VlX1u8CXgH/MdFxZ5CaRSDJm0s04JBLJxCOFQyKRZIwUDolEkjFSOCQSScZI4ZBIJBkjhUMikWSMFA6JRJIxUjgkEknGSOGQZB1FUdYoitKnKMpMRVGeUhRloupdJDlCZo5KcoKiKP8JVAH3CCEOTHQ8kuwihUOSExRFmY9VUr5UCNE90fFIsou8VJHkiuuwCg3/baIDkWQfKRySrKMoyreAXuAl4HpFUW6d4JAkWUZeqkgkkoyRMw6JRJIxUjgkEknGSOGQSCQZI4VDIpFkjBQOiUSSMVI4JBJJxkjhkEgkGSOFQyKRZMz/B9fnhtNWCm3WAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "train_data, = ax.plot(train_x.cpu().numpy(), train_y.cpu().numpy(), 'bo')\n",
    "\n",
    "model.eval()\n",
    "with torch.no_grad():\n",
    "    output = model.likelihood(model(train_x))\n",
    "    \n",
    "mean = output.mean\n",
    "lower, upper = output.confidence_region()\n",
    "line, = ax.plot(train_x.cpu().numpy(), mean.detach().cpu().numpy())\n",
    "ax.fill_between(train_x.cpu().numpy(), lower.detach().cpu().numpy(),\n",
    "                upper.detach().cpu().numpy(), color=line.get_color(), alpha=0.5)\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('y')\n",
    "ax.legend([train_data, line], ['Train data', 'Prediction'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next steps\n",
    "\n",
    "This was a pretty boring example, and it wasn't really all that different from GPyTorch's native SVGP implementation! The real power of the Pyro integration comes when we have additional latent variables to infer over. We will see an example of this in the [next example](./Clustered_Multitask_GP_Regression.ipynb), which learns a clustering over multiple time series using multitask GPs and Pyro."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
